Time Series Modeling - ARMA
自回归滑动平均模型ARMA(Autoregressive moving average model)
线形模型时间序列基础建模步骤
- 平稳性检验
首先对时间序列进行平稳性检验;若不平稳,需对时间序列做平稳化处理 平稳化方法:去除趋势(线性趋势,多项式趋势);差分;变换(对于方差变化的序列,经常采用log变化去除原来的指数趋势;更一般的,可以采用BOX-COX变换)
纯随机性检验(白噪声检验)
建模(定阶)
残差分析
模型优化(剔除不显著参数)
平稳性检验 — 单位根过程
建立ARMA模型的前提是时间序列是平稳的。
典型的非平稳时间序列是单位根(unitroot)非平稳时间序列。通常通过单位根检验检验时间序列是否平稳。
注:一些其他的平稳性检验方法有逆序数检验法,游程检验法,图检验方法。
典型的非平稳序列
随机游动 \[p_t=p_{t-1}+\epsilon_t\]
带漂移的随机游动 \[p_t=\mu+p_{t-1}+\epsilon_t\]
固定趋势模型 区别于上面两种,可以用回归/差分平稳。
ARIMA 通过一阶差分可以平稳。
单位根检验
单位根检验的理论基础:如果序列是平稳的,那么该序列的所有特征根应该在单位圆内。
DF检验
考虑如下的基础模型 \(x_t = \phi_0+\phi_1x_{t-1}+\xi_t, \xi_t\)为随机部分,且\(\xi_t\sim N(0,\sigma^2)\)
提出假设:原假设为序列非平稳\(H_0:\vert\phi_1\vert\ge1\),备择假设是序列平稳\(H_1:\vert\phi_1\vert\lt1\)
检验统计量: \[DF=\frac{\hat\phi_1-1}{SE(\hat\phi_1)}\] \(DF\)检验模型有三种: - 无漂移项自回归模型:\(x_t=\phi_1x_{t-1}+\xi_t\) - 有漂移项自回归模型:\(x_t = \phi_0+\phi_1x_{t-1}+\xi_t\) - 带趋势回归模型:\(x_t = \alpha + \beta t+\phi_1x_{t-1}+\xi_t\)
ADFTest
考虑如下的基础模型: \[X_t = c_t+\beta X_{t-1}+\sum_{j=1}^{p-1}\phi_j\Delta X_{t-j}+e_t\] 当\(\beta=1\)时,就是\(\Delta X_t\)的AR(p-1)模型,当\(\beta<1\)时,是\(X_t\)的AR(p)模型。 \(c_t\)是非随机的趋势部分,可以取0,或常数,或a+bt这样的非随机线性趋势。
检验假设: \(H_0:\beta=1\) vs. \(H_1:\beta<1\)
如果拒绝H0,就说明没有单位根,使用统计量\[ADF=\frac{\hat\beta-1}{SE(\hat\beta)}\] 当ADF统计量足够小的时候拒绝H0.
- fUnitRoots包的adfTest()函数可以执行单位根ADF检验。
- tseries包的adf.test()函数也可以执行单位根ADF检验。
注意,单位根DF检验和ADF检验都是在拒绝(显著)时否认有单位根, 不显著时承认有单位根。
PP检验:
ADF和DF检验要求随机扰动项独立同分布或者独立同正态分布,但是很多时间序列不满足这一点,故此提出\(PP\)检验,它限定随机序列可以具有暂时相关性和不同分布,它的原假设时间序列含有一个单位根,对于模型\(x_t = \phi x_{t-1}+\xi_t\). 提出原假设:\(H_0:\phi=1\)
ADF检验主要适用于方差齐性的场合,它对于异方差序列的平稳性检验结果不佳。
PP检验的适用于异方差场合的单位根检验,且PP检验的临界值和ADF检验相同,服从相应的ADF检验统计量的极限分布。
KPSS检验
KPSS检验:ADF检验和PP检验的原假设都是序列中有单位根,也就是认为序列是不平稳的,备择假设是序列是平稳的。在实证研究中,人们愿意多做差分(有单位根就需要做差分)而不愿意忽略单位根的存在,因为忽略单位根比过度差分带来的后果更为严重。其次以单位根作为原假设相对更加容易涉及检验统计量。故此以序列的平稳性为原假设的方法产生:KPSS方法。
单位根检验实例:
da <- read_table2("~/Downloads/ftsdata/q-gnp4710.txt", col_types=cols(
.default = col_double()))
gnp <- ts(log(da[["VALUE"]]),
start=c(1947, 1), frequency=4)
dgnp <- diff(gnp)
rm(da)
plot(dgnp, xlab="year", ylab="log(GNP)")##
## Call:
## ar(x = dgnp, method = "mle")
##
## Coefficients:
## 1 2 3 4 5 6 7 8
## 0.4318 0.1985 -0.1180 0.0189 -0.1607 0.0900 0.0615 -0.0814
## 9
## 0.1940
##
## Order selected 9 sigma^2 estimated as 8.918e-05
#AIC取9阶
#ADF的基础模型需要一个AR阶数,取p=9。
#用fUnitRoots::adfTest()对GNP的对数值进行ADF单位根检验:
library(fUnitRoots)
#type: c: with constant(intercept)
# nc: no constant
# ct:constant and trend
#这里gnp的均值非0,且存在trend,故type='ct'
fUnitRoots::adfTest(gnp, lags=9, type="ct") ##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 9
## STATISTIC:
## Dickey-Fuller: -0.0094
## P VALUE:
## 0.99
##
## Description:
## Tue Mar 9 00:19:35 2021 by user:
##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 9
## STATISTIC:
## DF: -0.0094
## P VALUE:
## t: 0.996
## n: 0.996
##
## Description:
## Tue Mar 9 00:19:35 2021 by user:
#结果值较大,说明不能拒绝零假设, 即对数GNP序列有单位根。序列非平稳
#另外两种单位根检验方法
#tseries::pp.test(gnp)
#tseries::kpss.test(gnp)
#尝试人为地拟合非随机线性增长趋势,检验残差是否有单位根
tmp.t <- c(time(gnp))
tmp.y <- residuals( lm(c(gnp) ~ tmp.t) )
fUnitRoots::adfTest(tmp.y, type="nc")##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## Dickey-Fuller: -0.0763
## P VALUE:
## 0.592
##
## Description:
## Tue Mar 9 00:19:35 2021 by user:
纯随机性检验(白噪声检验)
由\(Bartlett\)定理可以通过构造两种统计量:
\(Bartlett\)定理:如果一个序列是纯随机的,得到一个观察期数为\(n\)的观察序列\(x_t,t=1,2,...,n\),那么该序列的延迟非零期的样本自相关系数将近似服从均值为零,方差为序列观察期倒数的正态分布,即: \[ \hat\rho_k\sim N(0,\frac{1}{n}),\forall k\neq0 \] 由此提出原假设:
\(H_0:\rho_1=\rho_2=...=\rho_m=0\)(延迟期数小于等于\(m\)的序列值之间相互独立)
\(H_1:\)至少存在一个\(\rho_k\neq0\)
混成检验:
(1)\(Q\)统计量(适用于大样本): \[ Q_{BP}=n\sum_{k=1}^m\hat \rho_k^2\sim \chi^2(m) \] (2)\(LB\)统计量(适用于小样本) \[ Q_{LB} = n(n+2)\sum_{k=1}^m\frac{\hat\rho_k^2}{n-k}\sim\chi^2(m) \]
注: Box.test的lag怎么取?通常可以取ln(序列长度),min(20,序列长度−1)*
##
## Box-Ljung test
##
## data: gnp
## X-squared = 1213.4, df = 5, p-value < 2.2e-16
##
## Box-Ljung test
##
## data: gnp
## X-squared = 2110.2, df = 9, p-value < 2.2e-16
##
## Box-Ljung test
##
## data: gnp
## X-squared = 4228.2, df = 20, p-value < 2.2e-16
ARMA模型辨识
不同于AR,MA, ARMA模型并非根据acf,pacf定阶。
可以逐个从低阶模型尝试,p+q越小越好, 找到AIC最小的选择, 用精确最大似然或者条件最大似然方法估计参数。 对残差进行白噪声检验以验证模型是否充分。
R的forecast包提供了一个auto.arima()函数,可以自动进行模型选择。auto.arima()R中的函数使用了Hyndman-Khandakar算法的一种变体(Hyndman&Khandakar,2008),该算法结合了单位根检验,最小化AICc和MLE来获得ARIMA模型。auto.arima()为算法提供许多变体的参数。
例6.1考虑3M公司股票从1946年2月到2008年12月的月对数收益率, 共有755个观测。
d <- read_table2(
"~/Downloads/ftsdata/m-3m4608.txt",
col_types=cols(.default=col_double(),
date=col_date(format="%Y%m%d")))
mmm <- xts(log(1 + d[["rtn"]]), d$date) #计算月对数收益率
rm(d)
tclass(mmm) <- "yearmon"
ts.3m <- ts(coredata(mmm), start=c(1946,2), frequency=12)
head(ts.3m)## Feb Mar Apr May Jun
## 1946 -0.081125460 0.018421282 -0.105360516 0.190518702 0.005114897
## Jul
## 1946 0.073743834
ACF很接近于白噪声。PACF也比较接近于白噪声但是有比较多的超出界限的值, 尽管超出量不大。
## Series: ts.3m
## ARIMA(3,0,1)(1,0,1)[12] with non-zero mean
##
## Coefficients:
## ar1 ar2 ar3 ma1 sar1 sma1 mean
## 0.0453 -0.0285 -0.0837 -0.1124 0.5319 -0.4435 0.0103
## s.e. 0.3146 0.0417 0.0387 0.3147 0.2885 0.3049 0.0023
##
## sigma^2 estimated as 0.003998: log likelihood=1016.63
## AIC=-2017.25 AICc=-2017.06 BIC=-1980.24
auto.arima()函数选择了一个季节ARMA模型。
第一章作业:查找中国GDP季度数据,并尝试用ARIMA模型建模。
SARIMA
da <- read_excel('~/Desktop/应用时间序列实验课/实验课/数据/中国GDP数据.xlsx',sheet=2)
da <- ts(da$`季度GDP(亿元)`,frequency = 4,start=c(2008,1))
autoplot(da)自动拟合arima的阶数
## Series: da
## ARIMA(0,1,2)(1,1,0)[4]
##
## Coefficients:
## ma1 ma2 sar1
## -0.8613 0.2006 -0.4118
## s.e. 0.1209 0.1261 0.1140
##
## sigma^2 estimated as 33332300: log likelihood=-694.77
## AIC=1397.54 AICc=1398.17 BIC=1406.48
\[(1+0.4118B^4)(1-B)(1-B^4)y_t = (1-0.8613B+0.2006B^2)\epsilon_t\]
## Series: da
## ARIMA(0,1,2)(1,1,0)[4]
##
## Coefficients:
## ma1 ma2 sar1
## -0.8613 0.2006 -0.4118
## s.e. 0.1209 0.1261 0.1140
##
## sigma^2 estimated as 33332300: log likelihood=-694.77
## AIC=1397.54 AICc=1398.17 BIC=1406.48
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 580.3015 5452.413 2986.756 1.229971 3.998121 0.3203815 -0.01847175
##
## Box-Ljung test
##
## data: resm$residuals
## X-squared = 5.7441, df = 12, p-value = 0.9284
##
## Box-Ljung test
##
## data: resm$residuals
## X-squared = 11.725, df = 24, p-value = 0.9828
##
## Box-Ljung test
##
## data: resm$residuals^2
## X-squared = 13.657, df = 6, p-value = 0.03371
library(fGarch)
mod1 <- garchFit( ~ garch(1,1), data=resm$residuals, trace=FALSE)
#其中1表示均值方程是一个常数,
#输出结果中mu表示均值方程的均值,omega表示alpha0,alpha1为alpha1
resi <- residuals(mod1, standardize=TRUE)
checkresiduals(resi)##
## Box-Pierce test
##
## data: resi
## X-squared = 12.588, df = 12, p-value = 0.3997
##
## Box-Pierce test
##
## data: resi^2
## X-squared = 7.3936, df = 12, p-value = 0.8305
- 超前多步预报
tmp.y <- window(da, start=start(da), end=c(2024,4))
resm2 <- arima(tmp.y, order=c(0,1,2), seasonal=list(order=c(1,1,0), period=4))
#样本外预测
resm2 %>% forecast(h=12) %>% autoplot()mode2 <- garchFit( ~ garch(1,1), data=resm2$residuals, trace=FALSE)
pred1 <- predict(resm2, n.ahead=6)
mod <- predict(mode2,n.ahead=6)- 样本内静态预测(预测样本内最后6个观测)
horiz=6
dbnp.prev=predict(resm2,n.ahead=horiz)$pred
dbnp.se=predict(resm2,n.ahead=horiz)$se
bornesup=dbnp.prev+1.96*dbnp.se
borneinf=dbnp.prev-1.96*dbnp.se
library(ggplot2)
ts.plot(da)
lines(dbnp.prev,col="red")
lines(bornesup,col="orange")
lines(borneinf,col="orange")- 样本内动态预测(预测一个前进一个建模)
dbnp.prev.dyn=ts(start=length(da)-horiz+1,end=length(da))
dbnp.se.dyn=ts(start=length(da)-horiz+1,end=length(da))
for (i in (length(da)-horiz+1):length(da)){
dbnp.mod.temp=arima(window(c(da),end=i),order=c(0,1,2),seasonal=list(order=c(1,1,0),period=4),include.mean=F)
dbnp.prev.dyn[i-length(da)+horiz]=predict(dbnp.mod.temp,n.ahead=1)$pred
dbnp.se.dyn[i-length(da)+horiz]=predict(dbnp.mod.temp, n.ahead=1)$se
}
par(mfrow=c(1,1))
ts.plot(window(c(da),start=1,end=74),col="orange",main='Dynamic Prediction')
lines(dbnp.prev.dyn,col="red")
lines(dbnp.prev.dyn+1.96*dbnp.se.dyn,col="blue")
lines(dbnp.prev.dyn-1.96*dbnp.se.dyn,col="blue")## [1] 18377.35
ARIMA拟合
## Series: da
## ARIMA(4,1,1) with drift
##
## Coefficients:
## ar1 ar2 ar3 ar4 ma1 drift
## -0.0568 0.0093 0.0422 0.8713 -0.7753 1962.7763
## s.e. 0.0629 0.0571 0.0588 0.0538 0.0892 926.1329
##
## sigma^2 estimated as 38114765: log likelihood=-740.57
## AIC=1495.13 AICc=1496.86 BIC=1511.17
##
## Ljung-Box test
##
## data: Residuals from ARIMA(4,1,1) with drift
## Q* = 15.543, df = 3, p-value = 0.001407
##
## Model df: 6. Total lags used: 9
##
## Box-Pierce test
##
## data: arima.model2$residuals
## X-squared = 16.203, df = 12, p-value = 0.1821
##
## Box-Pierce test
##
## data: arima.model2$residuals^2
## X-squared = 17.662, df = 12, p-value = 0.1263
\[(1+0.0568B-0.0093B^2-0.0422B^3-0.8714B^4)(1-B)^1(y_t-1962t)=(1-0.7753B)e_t\]
非季节的ARIMA模型可以写为: \[(1-\phi_1B-...-\phi_pB^p)(1-B)^dy_t=c+(1+\theta_1B+...+\theta_qB)e_t\] 上式表达等价于: \[(1-\phi_1B-...-\phi_pB^p)(1-B)^d(y_t-\mu t^d/d!)=(1+\theta_1B+...+\theta_qB)e_t\] 其中B是滞后算子,\(c=\mu (1-\phi_1-...-\phi_p)\),及\(\mu\)是\((1-B)^dy_t\)的均值。 当d=0,且include.drift=T的时候,模型形式为: \[(1-\phi_1B-...-\phi_pB^p)(y_t-a-bt)=(1+\theta_1B+...+\theta_qB)e_t\] a:intercept;b:drift
- SARIMA 包含了季节项的滞后算子,如对于\(ARIMA(1,1,1)(1,1,1)_4\)模型, 模型可以写作: \[(1-\phi_1B)(1-\Phi_1B^4)(1-B)(1-B^4)y_t=(1+\theta_1B)(1+\Theta_1B^4)\epsilon_t\]
第二章课后练习查找格力电器(000651)日收盘价和收益率数据,分别检验其平稳性。
library(lubridate)
library(zoo)
library(xts)
library(tseries)
da <- read.csv('~/Desktop/000651.csv',fileEncoding = "GBK")
ts.stock <- zoo(da[,c('收盘价')],
as.POSIXct(da$日期))
plot(as.xts(ts.stock), type="l",
multi.panel=TRUE, theme="white",
major.ticks="years",
grid.ticks.on = "years")##
## Call:
## ar(x = as.xts(ts.stock))
##
## Coefficients:
## 1 2 3 4 5 6 7 8
## 0.5478 0.1864 0.0740 0.0377 -0.0087 0.0621 0.0122 0.0186
## 9 10
## 0.0267 0.0334
##
## Order selected 10 sigma^2 estimated as 11.55
##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 10
## STATISTIC:
## Dickey-Fuller: -1.3274
## P VALUE:
## 0.1933
##
## Description:
## Tue Mar 9 00:19:45 2021 by user:
##
## KPSS Test for Level Stationarity
##
## data: as.xts(ts.stock)
## KPSS Level = 14.588, Truncation lag parameter = 11, p-value = 0.01
##
## Phillips-Perron Unit Root Test
##
## data: as.xts(ts.stock)
## Dickey-Fuller Z(alpha) = -88.771, Truncation lag parameter = 11,
## p-value = 0.01
## alternative hypothesis: stationary
- 简单收益率
library(dplyr)
simple.return <- function(x){
x <- as.vector(x)
c(NA, diff(x) / x[1:(length(x)-1)])
}
d.geli <- data.frame(date=da$日期,收益率=simple.return(da$收盘价))
d.geli$收益率 <- as.numeric(d.geli$收益率)
d.geli <- d.geli[d.geli$收益率 != Inf & !is.na(d.geli$收益率),]
ts.plot(d.geli$收益率)ts.stock <- as.xts(zoo(d.geli$收益率,as.POSIXct(d.geli$date)))
fUnitRoots::adfTest(ts.stock,lags=10) # 平稳##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 10
## STATISTIC:
## Dickey-Fuller: -22.0079
## P VALUE:
## 0.01
##
## Description:
## Tue Mar 9 00:19:45 2021 by user:
##
## KPSS Test for Level Stationarity
##
## data: as.xts(ts.stock)
## KPSS Level = 0.30386, Truncation lag parameter = 10, p-value = 0.1
##
## Phillips-Perron Unit Root Test
##
## data: as.xts(ts.stock)
## Dickey-Fuller Z(alpha) = -5436.6, Truncation lag parameter = 10,
## p-value = 0.01
## alternative hypothesis: stationary
- 对数收益率
log.return <- function(x){
c(NA, diff(log(x)))
}
d.geli <- data.frame(date=da$日期,收益率=log.return(da$收盘价))
d.geli$收益率 <- as.numeric(d.geli$收益率)
summary(d.geli$收益率)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -Inf -0.01323 0.00000 NaN 0.01228 Inf 174
R Session Information
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.5
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] zh_CN.UTF-8/zh_CN.UTF-8/zh_CN.UTF-8/C/zh_CN.UTF-8/zh_CN.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] dplyr_1.0.2 tseries_0.10-47 lubridate_1.7.9
## [4] ggplot2_3.3.2 fGarch_3042.83.2 fUnitRoots_3042.79
## [7] fBasics_3042.89.1 timeSeries_3062.100 timeDate_3043.102
## [10] readr_1.4.0 readxl_1.3.1 xts_0.12.1
## [13] zoo_1.8-8 forecast_8.13
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.1.0 xfun_0.18 purrr_0.3.4 urca_1.3-0
## [5] lattice_0.20-41 colorspace_1.4-1 vctrs_0.3.4 generics_0.0.2
## [9] htmltools_0.5.1.1 yaml_2.2.1 rlang_0.4.10 pillar_1.4.6
## [13] withr_2.4.1 glue_1.4.2 TTR_0.24.2 lifecycle_1.0.0
## [17] quantmod_0.4.17 stringr_1.4.0 munsell_0.5.0 gtable_0.3.0
## [21] cellranger_1.1.0 evaluate_0.14 labeling_0.4.2 knitr_1.30
## [25] rmdformats_0.3.7 lmtest_0.9-38 parallel_4.0.3 curl_4.3
## [29] Rcpp_1.0.6 scales_1.1.1 farver_2.0.3 fracdiff_1.5-1
## [33] hms_0.5.3 digest_0.6.27 stringi_1.5.3 bookdown_0.21
## [37] grid_4.0.3 quadprog_1.5-8 tools_4.0.3 magrittr_2.0.1
## [41] tibble_3.0.4 crayon_1.4.1 spatial_7.3-12 pkgconfig_2.0.3
## [45] ellipsis_0.3.1 rmarkdown_2.5 R6_2.5.0 nnet_7.3-14
## [49] nlme_3.1-149 compiler_4.0.3